// -------------------------------------------------------------------------
// -----                 EiccFdtIdealRecoTask header file             -----
// -----                  Created 2021/04/30  by Y. Liang             -----
// -------------------------------------------------------------------------


/** EiccFdtIdealRecoTask.h
 *@author Y. Liang <liangyt@impcas.ac.cn>
 **
 ** The ideal tracking produces tracks of type FdtPoint by using MC information
 ** directly, but with the correct data interface.
 **/


#ifndef EICCFDTIDEALRECOTASK_H
#define EICCFDTIDEALRECOTASK_H

// Root includes
#include "TVector3.h"
#include "TRandom.h"
#include "TString.h"
#include "TGeoManager.h"

// framework includes
#include "FairTask.h"

// EiccFdt includes
#include "EiccFdtPoint.h"
#include "EiccFdtHit.h"
#include "EiccMCTrack.h"
//#include "EiccGeoHandling.h"


//Genfit2
#include <AbsFinitePlane.h>
#include <AbsFitterInfo.h>
#include <AbsMeasurement.h>
#include <AbsTrackRep.h>
#include <ConstField.h>
#include <DetPlane.h>
#include <Exception.h>
#include <FieldManager.h>
#include <KalmanFittedStateOnPlane.h>
#include <AbsKalmanFitter.h>
#include <KalmanFitter.h>
#include <KalmanFitterRefTrack.h>
#include <KalmanFitterInfo.h>
#include <KalmanFitStatus.h>
#include <DAF.h>
#include <GFGbl.h>
#include <MeasuredStateOnPlane.h>
#include <MeasurementOnPlane.h>
#include <FullMeasurement.h>
#include <PlanarMeasurement.h>
#include <ProlateSpacepointMeasurement.h>
#include <RectangularFinitePlane.h>
#include <ReferenceStateOnPlane.h>
#include <SharedPlanePtr.h>
#include <SpacepointMeasurement.h>
#include <StateOnPlane.h>
#include <Tools.h>
#include <TrackCand.h>
#include <TrackCandHit.h>
#include <Track.h>
#include <TrackPoint.h>
#include <WireMeasurement.h>
#include <WirePointMeasurement.h>

#include <MaterialEffects.h>
#include <RKTools.h>
#include <RKTrackRep.h>
#include <StepLimits.h>
#include <TGeoMaterialInterface.h>

#include <EventDisplay.h>

#include <HelixTrackModel.h>
#include <MeasurementCreator.h>


class EiccFdtHit;
class TClonesArray;

class EiccFdtIdealRecoTask : public FairTask
{
 public:

  /** Default constructor **/
  EiccFdtIdealRecoTask();

  EiccFdtIdealRecoTask(Double_t sx, Double_t sy, Double_t sz, bool debug, TString foutname, TString fouteffname, TString fbtype);


  /** Destructor **/
  virtual ~EiccFdtIdealRecoTask();

  EiccFdtIdealRecoTask(const EiccFdtIdealRecoTask& ) = delete;
  EiccFdtIdealRecoTask& operator=(const EiccFdtIdealRecoTask& ) = delete;

  /** Virtual method Init **/
  virtual void SetParContainers();
  virtual InitStatus Init();


  /** Virtual method Exec **/
  virtual void Exec(Option_t* opt);

  virtual void FinishTask();

  void SetBranchName(TString name){fBranchName = name;};
  void SetOutputName(TString name){fOutname = name;};
  void SetOutputEffName(TString name){fOutEffname = name;};
  void SetBtype(TString name){fBtype = name;};
  void SetPdg(int input_pdg){pdg = input_pdg;};
  void SetDebug(bool debugkey){fDebug = debugkey;};
  void SetEvtDis(bool evtdis){DisplayEvent = evtdis;};
  void SetPersistance(Bool_t p = kTRUE){fPersistance=p;};
  Bool_t GetPersistance() {return fPersistance;};
  
private:
  Bool_t fPersistance; // switch to turn on/off storing the arrays to a file
  
  TString fBranchName;
  //EiccGeoHandling* fGeoH;
  /** Input array of Points **/
  TClonesArray* fPointArray;
  TClonesArray* fMctruthArray;

  /** Output array of Hits **/
  TClonesArray* fHitOutputArray;
  /** Properties **/
  Double_t fSigmaX;       // Variance vector of sensor smearing
  Double_t fSigmaY;       // Variance vector of sensor smearing
  Double_t fSigmaZ;       // Variance vector of sensor smearing
  EiccFdtPoint* fCurrentEiccFdtPoint;
  EiccMCTrack* fCurrentEiccMCTrack;
  TGeoHMatrix* fCurrentTransMat;
  TMatrixT<Double_t> fHitCovMatrix;
  bool fDebug;
  /** smearing and geometry access **/
  void InitTransMat();
  void smear(TVector3& pos);
  void smearLocal(TVector3& pos);
  void CalcDetPlane(TVector3& o,TVector3& u,TVector3& v);

  int randomPdg();
  int randomSign();
  void smearpixel(double pixelsize, int det_id, double input_X, double input_Y, double input_Z, double &output_X, double &output_Y, double &output_Z);


  void Register();

  void Reset();

  void ProduceHits();



  //GenFit variables
  const int nMaxHit = 100;
  const double maxmomentum = 5.0;     // GeV
  double pixel_size = 0.001;   // cm; size of the pixel 10 micron
  const double resolution = 0.001;   // cm; resolution of generated measurements
  const unsigned int nMeasurements = 6;
  double BField = 15.;       // kGauss
  int pdg = -11;               // particle pdg code
  const bool matFX = true;         // include material effects; can only be disabled for RKTrackRep!
  double momentum = 0.1;     // GeV
  const bool smearOn = true;

  const double thetaDetPlane = 90;         // degree
  const double phiDetPlane = 0;         // degree
  const double resolutionWire = 5*resolution;   // cm; resolution of generated measurements
  TVector3 wireDir;
  const double skewAngle = 5 ;
  const bool useSkew = true;
  const int nSuperLayer = 10;
  const double minDrift = 0.;
  const double maxDrift = 2;
  const bool idealLRResolution = false; // resolve the l/r ambiguities of the wire measurements

  const double outlierProb = -0.1;
  const double outlierRange = 2;

  const double hitSwitchProb = -0.1; // probability to give hits to fit in wrong order (flip two hits)

  const int splitTrack = -5; //nMeasurements/2; // for track merging testing.
  const bool fullMeasurement = false;

  const genfit::eFitterType fitterId = genfit::RefKalman;
  const genfit::eMultipleMeasurementHandling mmHandling = genfit::unweightedClosestToPrediction;

  const int nIter = 20; // max number of iterations
  const double dPVal = 1.E-3; // convergence criterion

  const bool resort = false;
  const bool prefit = false; // make a simple Kalman iteration before the actual fit
  const bool refit  = false; // if fit did not converge, try to fit again

  const bool twoReps = false; // test if everything works with more than one rep in the tracks

  const bool checkPruning = true; // test pruning

  const bool smearPosMom = true;     // init the Reps with smeared pos and mom

  const double posSmear = 10*resolution;     // cm
  const double momSmear = 5. /180.*TMath::Pi();     // rad
  const double momMagSmear = 0.2;   // relative
  const double zSmearFac = 2;

  //const bool debug = false;
  const bool onlyDisplayFailed = false; // only load non-converged tracks into the display
  //const bool DisplayEvent = false; // only load non-converged tracks into the display
  bool DisplayEvent = false; // only load non-converged tracks into the display

  const double charge = TDatabasePDG::Instance()->GetParticle(pdg)->Charge()/(3.);
/*
  double X[nMaxHit],Y[nMaxHit],Z[nMaxHit];
  double PX[nMaxHit],PY[nMaxHit],PZ[nMaxHit];
  int DetID[nMaxHit];
*/
  double X[100],Y[100],Z[100];
  double PX[100],PY[100],PZ[100];
  int DetID[100];
  
  double TruthX[100],TruthY[100],TruthZ[100];
  double TruthPX[100],TruthPY[100],TruthPZ[100];
  int TruthMotherID[100];

  std::vector<genfit::eMeasurementType> measurementTypes;

  genfit::AbsKalmanFitter* fitter = 0;
  genfit::MeasurementCreator measurementCreator;

  genfit::EventDisplay* display;

  genfit::Track* fitTrack;
  genfit::Track* secondTrack;

  TString fOutname;
  TString fOutEffname;
  TString fBtype;
  TFile *f_data_out;
  TTree *t_out;
  
  double px_truth, py_truth, pz_truth, pt_truth, p_truth, theta_truth, phi_truth, eta_truth;
  double px_rec, py_rec, pz_rec, pt_rec, p_rec, theta_rec, phi_rec, eta_rec;
  double mom_res, mom_res_old, dca_rphi, dca_z;
  // double x_truth[100], y_truth[100], z_truth[100], r_truth[100], phi_truth[100], det_id[100], Evt_id[100];
  // double x_rec[100], y_rec[100], z_rec[100], phi_rec[100];
  double chi2;
  double abs_dca_rphi, dca_x, dca_y;
  int NHit;
  int fEventNo;

  std::ofstream fout;
  std::ofstream fouteff;

  ClassDef(EiccFdtIdealRecoTask,1);

};

#endif
